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SIMULATION  OUTPUT  ANALYSIS  FOR 
GENERAL  STATE  SPACE  MARKOV  CHAINS 

by 

Peter  W.  Glynn  and  Donald  L.  Iglehart 
Stanford  University 

1.  Introduction 

The  statistical  analysis  of  simulation  output  has  been  the  primary  focus 
of  recent  research  in  simulation  methodology.  Methods  have  been  developed 
which  permit  the  simulator  to  construct  confidence  intervals  for  steady-state 
characteristics  of  the  system  being  simulated.  The  principal  methods  in 
current  use  are  autoregressive  modeling,  batch  means,  regenerative,  and 
replication.  With  the  exception  of  the  replication  method,  all  methods  are 
based  on  just  one  simulation  run.  These  methods  for  constructing  confidence 
intervals  are  all  based  on  central  limit  theorems  for  the  underlying  stochastic 
processes  being  simulated.  Thus  all  methods  are  only  valid  asymptotically 
for  long  simulation  runs. 

In  this  paper  we  discuss  three  new  methods  for  analyzing  simulation 
output.  All  three  methods  are  aimed  at  analyzing  the  simulation  output  of 
general  state  space  Markov  chains.  This  class  of  processes  encompasses  the 
embedded  jump  chain  generated  by  a  generalized  semi-Markov  process  (GSMP) . 
GSMP's  are  important  for  simulation  since  they  may  be  used  to  model  a  general 
discrete  event  simulation.  GSMP's  have  been  discussed  in  recent  papers  Dy 
FOSSETT  (1979),  HORDIJK  and  SCHASSBERGER  (1981),  and  WHITT  (1980). 

The  first  method,  called  the  extended  regenerative  method  (ERM) ,  is  based 
on  some  recent  work  on  general  state  space  Markov  chains  by  ATHREYA  and  NEY 
(1978)  and  NUMMELIN  (1978).  This  method  involves  a  construction  which  creates 


regeneration  points  for  Harkov  chains  which  do  not  hit  a  single  point  infinitely 
often.  While  this  idea  is  very  attractive  in  principle,  there  are  a  number 
of  practical  considerations  which  limit  its  application.  However,  the 
method  can  be  used  to  increase  the  rate  of  regeneration  points  when  using 
the  standard  regenerative  method.  For  more  details  on  chis  method  and  some 
related  results  see  GLYNN  (1981). 

The  second  method,  called  the  random  blocking  method  (RBM) ,  is  based  on 
blocks  of  the  process  which  begin  when  the  process  enters  a  given  set  in 
the  state  space.  This  method  is  reminiscent  of  the  regenerative  method  except 
the  blocks  created  here  are  not  independent  and  identically  distributed. 

Details  on  this  second  method  can  be  found  in  a  forthcoming  paper  by  GLYNN 
(1981). 

The  last  method  is  a  variation  of  the  method  of  autoregressive  modeling. 
This  method,  the  multivariate  autoregressive  method  (MARM) ,  fits  a  multi¬ 
variate  autoregressive  model  to  the  simulation  output  data.  The  model  fitting 
is  done  automatically  based  on  Akaike’s  AlC-criterion;  see  AKAIKE  (1976)  for 
a  full  discussion  of  these  criterion.  A  forthcoming  paper  by  JOW  (1981)  will 
develop  this  method  and  other  related  methods  for  simulation  applications. 

This  paper  is  organized  as  follows.  Section  2  is  denoted  to  a  discussion 
of  GSMP's  and  their  relation  to  simulation.  The  ERM  is  covered  in  Section  3 


and  RBM  in  Section  4. 
illustrative  example 
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2.  Generalized  Semi-Markov  Processes 

In  a  discrete-event  simulation  a  finite  number  of  events  occurring  at 

random  times  cause  changes  in  the  state  of  the  system  being  simulated.  The 

number  of  events  active  at  any  given  time  is  a  function  of  the  state  of  the 

system.  This  type  of  simulation  is  well  modeled  by  the  generalized  semi- 

Markov  process  (GSMP)  which  we  now  describe. 

Let  S  be  the  finite  (or  countable)  set  of  states  which  describes  the 

GSMP  at  the  successive  transition  epochs  and  G  the  finite  number  of  events 

which  can  cause  a  transition.  If  G  «  {e, ,  e„,  ...»  e  },  then  we  let 

1  Z  m 

G(s)  =  (e,(s),  ...,  e  (s)}  denote  the  subset  of  events  active  when  the  GSMP 
l  n 

s 

is  in  state  s.  For  each  event  active  in  state  s,  associate  a  clock  which 
records  the  time  until  that  event  would  trigger  a  state  change.  If  in 
state  s,  the  clock  (associated  with  an  event  in  G(s))  with  the  minimal 
reading  triggers  the  next  state  change  when  it  runs  down  to  zero.  Let  C(s) 
denote  the  possible  clock  readings  when  the  GSMP  is  in  state  s: 

C(s)  =  {c  €  1R™  :  c^  >  0  iff  e^  €  G(s);  c^  ^  c^  for  i  t  j}  , 

where  ]R™  is  the  Cartesian  product  of  m  copies  of  [0,°°).  Next  define  the 
space  2  by 


2  =  U  ({s>  x  C(s)) 
s  €  S 


and  the  process  (S,C)  =  { (S  ,  C  )  :  n  >  0}  which  lives  on  2  and  represents 
the  state  values  and  clock  readings  at  the  successive  transition  epochs.  The 
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process  (£,£)  is  a  general  state  space  Markov  chain  whose  transition  kernel 
will  be  defined  in  Section  3.  Finally,  the  GSMP  is  a  piece-wise  constant 
process,  {X  :  t  _>  0),  constructed  from  the  embedded  jump  process  (S,C)  in 
the  usual  manner. 


3.  The  Extended  Regenerative  Method  for  GSSMC's 

We  start  by  formalizing  the  notion  of  a  general  state  space  Markov  chain 
(GSSMC).  Let  E  be  a  complete,  separable  metric  space  with  E  its  associated 
Borel  field.  A  function  P  :  E  x  E  [0,1]  is  called  a  probability  transition 
kernel  if : 

i)  P(x,*)  is  a  probability  on  (E,E)  for  each  x  in  E, 
ii)  P(-,B)  is  E-measurable  for  all  B  €  E. 

One  should  think  of  P(x,B)  as  representing  the  one-step  transition  probability 
of  the  chain  passing  from  x  into  B.  The  analogous  n-step  transition  prob¬ 
abilities  are  then  given  through  the  Chapman-Kolmogorov  equations,  namely 

Pn+1(x,B)  =  /  Pn(y,B)  P (x ,dy) 

E 

where,  of  course,  P^(x,B)  =  5  (B)  (6  (B)  is  1  or  0  depending  on  whether 

X  X 

or  not  x  €  B) .  Given  a  kernel  P  and  an  initial  probability  p  on 

(E,E),  one  can  construct  a  measure  P  on  (2,F)  =  (E  x  E  x  • • • ,  E  x  E  x  • • • ) 

such  that 

pBfxo  *  V  *1  <  "i-  {  V 

=  /  P.(dxn)  /  P(xn,  dx  )  /  P(x  .  ,  dx  ) 

Bn  °  B.  °  1  B  n_1  n 

0  1  n 

k 


where  X^(cu)  =  and  (co^,  co^,  ...)  is  a  typical  element  of  S2.  The 
above  construction  yields  a  process  {X^  :  n  J>  1}  that  is  endowed  with  the 
Markov  property 


Wi e  b|xo . v 


P  {X  , ,  <E  B|x  } 
p.  n+1  n 


and  is  referred  to  as  the  GSSMC  associated  with  kernel  P  and  initial 
distribution  p.. 


EXAMPLE  3.1.  Let  E  =  {0,  1,  ...}  and  put  P(i,B)  =  J,.  p  . ,  where  P  is 

J  ij 

a  stochastic  matrix  on  E.  This  gives  the  classical  countable  state  Markov 
chain. 


EXAMPLE  3.2.  Let  (W^  :  n  J>  0}  be  the  waiting  time  process  of  a  GI/G/1 
queue.  Suppose  that  u  and  v  are  the  interarrival  and  service  r.v.'s 
respectively,  and  that  v-u  has  a  density  f  (x) .  Then  E  =  and 

x 

P(x,B)  =*  /  f  (y-x)dy  +  5  (B)  /  f  (y)dy  . 

B  0 

EXAMPLE  3.3.  Consider  a  model  that  imitates  the  dynamics  of  a  lake.  Let  S^, 
and  X^  represent  the  volume  of  water  stored  in  the  lake,  the  inflow  of 
water,  and  outflow  of  water  respectively,  at  time  i.  Then,  the  mass-balance 
equation 


(3.4) 


S 


i-1 


4 
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holds.  If  one  now  assumes  that  output  increases  linearly  with  storage  through 
the  relation  =  aS^  (0  <  a  <  1),  then  (3.4)  takes  the  form 


^i-l  +  £i 


where  p  =  1/ (1+a)  and  =  (aZ^)/(l+a).  Finally,  the  further  assumption 
that  {e^}  is  i-i-d.  with 

P{e  €  B}  =  P60(B)  +  (1-p)  /  f (y)dy 

B 

yields  a  Markov  chain  model  for  the  outflows  {X^},  where  E  =  [0,°°)  and 


P(x,B)  =  p6  (B)  +  (1-p)  /  f (y-px)dy  . 

B 

EXAMPLE  3.5.  The  embedded  jump  process  (S,C)  introduced  in  Section  2  is  a 
GSSMC  on  state  space  Z  and  with  kernel 

P((s,c) ,  B) 

*  p(s' ;  s,  i*)  17  /  f(y;  s',i,s,i*)dy  n  6  .(B) 

i€Ns,  Bi  j  €0g,  Ci  J 

where  B  is  that  subset  of  Z  corresponding  to  the  GSMP  entering  state  s' 
with  the  i'th  clock  of  s'  set  to  a  value  in  B.. 

i 

The  numbers  p(s' ;  s,i*)  govern  state  transitions  of  the  GSMP  and  represent 

the  probability  of  a  jump  to  s'  from  s,  given  that  clock  i*  initiated 

the  jump.  The  rest  of  the  kernel  governs  clock  readings.  Those  clocks 

i  €  N  ,  (now  clocks)  are  set  stochastically  according  to  a  density 
s 
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f (y ;  s',i,s,i*),  whereas  those  clocks  j  €  0gl  (old  clocks)  are  inherited  from 

the  previous  state  and  so  must  be  set  deterministically  at  the  previous  value 

c*. 

J 

Note  that  by  setting  x  =  c,(i*(C,))  +  •••  +  c  (i*(C  ))  (c.(k)  is  the 

'  nil  n  n  j 

value  of  the  k'th  clock  at  the  j’th  jump  of  X£) .  We  can  retrieve 

{X_  :  t  >  0>,  the  GSMF,  from  { (S  ,C  )  :  n  >  0}  via 
t  —  ’  n  n  — 

oo 

xc  ■  l  6c(1V  Vl»  Sk  • 

k*0 

Let  us  turn  now  to  a  recurrence  condition  for  GSSMC's  first  formulated 

by  ATHREYA  and  NEY  (1978)  and  NUMMELIN  (1978).  We  say  that  (Xn  :  n  >  0} 

is  (A, B,\,cp ,k)  recurrent  if  there  exist  A,  B  €  E,  a  positive  number  X, 

an  integer  k,  and  a  probability  cp  on  B,  such  that: 

i)  p{£n=0  6X  (A)  =  +O0lxo  =  x)  *  1  for  all  x  in  E, 
n 

ii)  P  (x,E)  2  ^P(E)  for  each  x  in  A  and  measurable  subset  E  of  B. 

In  our  setting,  this  recurrence  notion  is  in  fact  equivalent  to  one  first 
proposed  by  HARRIS  (1956) . 

To  gain  some  appreciation  for  the  significance  of  these  conditions,  we 
observe  that  in  the  (A,B,X,<p,l)  case,  we  can  decompose  P  over  A  as 

P(x,E)  =  Vp(E)  +  (1-X)  Q(x,E) 

where  Q(x,*)  is  a  probability  on  (E,E).  The  key  idea  of  Athreya,  Ney  and 
Nummelin  was  to  exploit  this  decomposition  via  the  following  embedding. 

Let  E'  =  E  *  {0,1}  and  E'  be  the  associated  product  a-field.  We 
extend  P  to  a  kernel  P'  on  (E',E')  by  setting 
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1 


(  (1-X)  P (x ,E) ,  x<A 
P’ ((x,S),  E  x  {0})  *  { 

{  (1-X)  Q(x,E),  x  €  A 


(  XP(x,E) ,  x  f  A 

P'((x,6),  E  x  {1})  =  J 

'  Xp  (E  D  B) ,  X  €  A 


For  a  probability  p  on  (E ,  E),  define  p*  on  (E',E’)  by 


p' (E  x  {i»  =  60(i)  (1-X)  n(E)  +  61(i)X  n(E) 


Then,  it  can  be  readily  verified  that  the  Markov  chain  X'  *  (X  ,  5  )  on 

n  n  n 

(E',E')  associated  with  P* ,  p'  has  the  property  that  the  coordinate  process 
{X^  :  n  ^  0}  is  the  original  GSSMC  on  (E,E)  associated  with  P,  p. 

The  importance  of  this  embedding  is  that  it  furnishes  one  with  a 
sequence  of  regeneration  times  T^  defined  in  terms  of  X^  via 


T,  *  inf {n  >  1  :  X  .  €  A,  6  *  1}  , 

1  —  n-1  n 

T.  *  inf{n  >  T.  ,  :  X  ..  6  A,  6  =1}  ,  k>l. 
k  k-1  n-1  n  ’ 


The  regenerative  character  of  Xr  guarantees  "nice"  ergodic  behavior  for 

n 

X  .  For  example,  if  we  assume  that  E'(T.-T..)  <  4",  then  the  classical 
n  i  i 

regenerative  theory  shows  that  (e.g.,  SMITH  (1955)) 


\  j,  £<\>  **■  ■  V® 


a.  s. 


where 
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1 


T  -1 
2 


tt(B)  =  E'C  l  6X  (B)]  /E'(T2-T1) 


j=T,  J 


Estimation  of  r  is  a  common  goal  of  simulators.  The  above  discussion 
suggests  an  "extended  regenerative  method"  (ERM)  for  producing  confidence 
intervals  for  r  based  on  the  sequence. 

1.  Generate  the  sequence  { (X  ,  6  )  :  n  >  0}. 

n  n  — 

Tk+1 

2.  Let  Yk  =  l  f(Xj),  ak  =  ?k+1  -  Tk  and  form 

j=Tk 

?(n)  =  C  l  \)/C  l 

k=l  k=l 

Sll(n)  =  n-1  Yk  -  n(n-l)  YlP 

S12(n)  =  Yk  ak  "  n(n-l)  ^  YlP  ^  ak^ 

s~-  (n)  =  — -  £  a, - — r-r  ^  \  a  ") 

22  n-1  k“^  k  n(n-l)  ^  kf^  lr 


s2(n)  =  s^Cn)  -  2r(n)  s12(n)  +  f2(n)  s22(n)  . 


3.  To  form  a  100(1-6)%  confidence  interval  for  r  choose  z  so  that 

6 

$(zc)  =  1  -  6/2,  where  $>  is  the  standard  normal  distribution  function, 
o 


Then 


z  s (n)  n1/2  „  z  s (n)  n1/2  n 

r(n)  _  ,  r(n)  +  0 


n 

l  \ 

k=l  K 

is  the  desired  confidence  interval. 


n 

l  \ 

k=l  K 
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The  above  steps  can  all  be  justified  under  the  assumption  that 
2 

0  <  E(Y^  -  r a^)  <  +°°.  Incidentally,  although  Step  1  appears  to  necessitate 

the  ability  to  generate  deviates  distributed  according  to  Q (x , • ) »  this 
difficulty  can  be  avoided  by  using  an  acceptance-rejection  technique  ([8],  p.  57). 


EXAMPLE  3.1  (continued).  Suppose  that  X  =  {X^  :  n  _>  0}  is  an  irreducible, 
recurrent  Markov  chain  with  transition  matrix  P  =  {p^}.  Then,  X  is 
({0},  E,  1,  Pq  ,1)  recurrent  and  T^  =  N^+1,  where  is  the  k'th  hitting 

time  of  0.  The  ERM  reduces  here  to  the  classical  regenerative  technique 


(modulo  a  shift  to  the  right). 


EXAMPLE  3.2  (continued).  Assuming  that  Ev  jC  Eu,  {W^  :  n  0}  is  ({0},  E,  1, 
P (0, • ) ,  1)  recurrent,  and  again  T^  =  N^+1,  when  is  the  k'th  hitting 

time  of  0,  reducing  our  estimation  technique  to  the  usual  one  (again, 
modulo  a  shift). 


EXAMPLE  3.3  (continued).  Assuming  that  Es^  <  +°°  and  that  f  is  positive 

and  continuous  over  [0,  +»),  {X  :  n  >  0)  is  ( [ 0 , b] ,  E,  X,  cp ,  1)  recurrent 

n  — 

where 

My)  3  (1-p)  min  f(y-qx) 

0  _<  x  <  b 

CO 

x  =  /  t(y)dy;  <p(y)  =t(y)A  . 
o 

Here  the  T  ' s  form  a  subsequence  of  the  hitting  times  of  [0,b]  —  this 
reflects  the  fact  that  X^  returns  to  no  point  infinitely  often.  The  ERM 
is  applicable  here,  whereas  the  classical  regenerative  method  is  not. 


i 


EXAMPLE  3.5  (continued).  We  say  that  the  GSMP  has  a  single  set  if  there 
exists  a  state  s'  with  which  is  associated  only  one  clock.  Then,  given 
that  the  GSMP  hits  {s’}  infinitely  often,  (S,C)  is  (A,  Z,  1,  P((s',c),  •),  1) 
recurrent  where  A  =  {(s,c)  €  Z  :  s  =  s'}.  Here,  =  N^+l  where  the  N^'s 
are  consecutive  hitting  times  of  A. 

It  should  be  remarked  that  in  the  context  of  a  GSMP,  a  simulator  is 
often  interested  in  the  continuous-time  quantity 

t 

"r  =  lira  /  f  (X  )ds/t 
t  ■*■  00  0  s 


This  can  be  estimated  by  the  ERM,  provided  that  the  definitions  of  and 

a.  are  modified  to 
k 


£<V  V“<cj” 


rw 


l  c  (i*(C  )) 


The  above  discussion  focussed  on  (A,  B,  X,  cp ,  1)  GSSMC’s.  A  fundamental 
difficulty  arises  in  the  (A,  B,  X,  cp,  k)  case,  a  difficulty  that  can  be 
partially  circumvented  via  an  embedding  that  endows  the  chain  with  an 
environment  that  is  "loosely  regenerative”  in  the  sense  of  SMITH  (1955).  An 
estimation  procedure  can  then  be’  developed  that  bears  close  resemblance 
to  the  "extended  regenerative  method"  outlined  above  (see  [8]). 

Before  leaving  this  topic,  it  should  be  mentioned  that  a  number  of 
difficulties  remain  in  terms  of  implementation  of  the  ERM.  The  most  fundamental 
problem  is  that  the  decomposition  of  P  over  A  requires  an  explicit  form 
for  the  kernel  over  that  set.  The  "event-scheduling"  approach  normally  used 


11 


by  simulators  to  generate  sample  paths  does  not  require  an  explicit 
representation  for  the  kernel,  and  so  a  simulator  is  left  with  the  burden¬ 
some  task  of  calculating  such  a  form.  Additional  difficulties  arise  in 
determining  appropriate  X,  (p,  and  A,  although  here  simple  numerical 
techniques  would  be  applicable  ([8],  p.  52). 


4.  Random  Blocking  Method  for  GSSMC's 

Our  major  incentive  in  studying  GSSMC's  here  has  been  in  terms  of 
application  to  GSMP's.  It  turns  out  that  GSMP's  possessing  no  single  set 
cannot  be  (A,  B,  X,  <p,  1)  recurrent  ([8],  p.  60),  and  hence  the  ERM 
discussed  in  the  previous  section  does  not  apply.  This  factor  together 
with  the  ERM-related  difficulties  already  mentioned,  motivates  development 
of  other  methods. 

Suppose  that  {Xfl  :  n  _>  0}  is  an  (A,  B,  X,  u> ,  k)  GSSMC  with  invariant 
probability  tt  .  This  will  in  fact  be  the  case  for  the  embedded  jump  process 
(S,C),  corresponding  to  a  GSMP  {Xt  :  t  >  0},  provided  that  (see  [9]) 

i)  the  state  space  S  is  finite 

ii)  the  "road  map"  p(s',  s,  i*)  satisfies  a  natural  irreducibility 
condition  (see  [7],  p.  16) 

iii)  the  densities  f(’5  s',  i,  s,  i*)  are  positive  and  continuous 
on  (0,°°),  with  finite  mean. 

For  such  an  (A,  B,  X,  (p,  k)  recurrent  {X  :  n  >  0}  one  can  show  that 

n  — 

if  Ejf  (X)  |  <  +»,  then 

j0  £<V  *  r  ■  V® 
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a.  s. 


for  any  p.  on  (E,E).  A  simulator  commonly  wishes  to  obtain  an  estimate 
for  r,  together  with  an  associated  confidence  interval. 

Let  T^,  ^2*  be  the  consecutive  hitting  times  of  the  set  A.  The 

"random  blocking  method"  (RBM)  hinges  on  the  observation,  due  essentially  to 
OREY  (1959),  that 


i  }  )  l  O  I  •  *  •  9  ) 

k  *k  k+1 


is  a  Doeblin  recurrent  Markov  chain  with  a  single  ergodic  set  (see  DOOB 
(1953)  for  definitions  and  results).  Under  the  assumption  that  is 

aperiodic  (this  will  generally  be  the  case  for  GSMP's),  functions  of 
will  enjoy  a  central  limit  theorem  with  a  variance  constant  of  the  form 


(4.1) 


ct2  =  a2(Zn)  +  2  1  cov(ZQ,  Zk) 


k=l 


By  truncating  the  infinite  sura  and  estimating  the  finite  number  of  remaining 
terms  in  the  series  by  the  standard  sample  moments,  we  obtain  the  RBM. 

1.  Choose  a  truncation  number  ra  (the  number  of  covariance  terms  of 

(4.1)  to  be  retained). 

T  -1 

2.  Put  f<V-  \  *  Vl  -  Tt  a"d£°™ 


r(n)  -  C  I  O/C  I  \) 
k=l  k-1 

n-£ 


Co£(tl)  =  n -l  Yk  Yk U~  nV  \  Yk+1 


r(n) 


n-l 


*2 .  »  n-£ 

+  — —  l  a  a 


n-6  ,  ^  Yk  \+l  '  n-t  ,  L .  "k  “k+f 

k=l  k=l 
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/Eli  1  f(X  )|2+5>ir(dx)  <  +~  , 

A  X  j=l  2 

for  some  5  >  0  ([9]).  The  GSlfP's,  the  RBM  can  be  modified  in  the  same  way 

as  the  ERM  so  as  to  provide  confidence  intervals  for  continuous-time  quantities 
of  the  form  (3.5). 

Observe  that  when  m  =  0,  this  technique  reduces  to  the  "approximate 
regenerative  method"  (ARM)  of  CRANE  and  IGLEHART  (1975).  Hence  one  can  think 
of  the  RBM  as  a  second-order  refinement  of  the  ARM.  It  should  be  noted, 
however,  that  the  RBM  involves  the  undesirable  element  of  having  to  make  an 
a  priori  judgment  as  to  an  appropriate  value  of  the  truncation  index  m. 

This  is  a  problem  to  which  we  intend  to  devote  more  attention. 
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Multivariate  Autoregressive  Method 

The  last  method  we  shall  discuss  is  the  multivariate  autoregressive 

method  (MARM) .  Let  {X  :  n  >  0}  be  a  vector-valued  strictly  stationary 

~n  — 

process  with  mean  vector  =  E{Xq}  and  covariance  function  £(h) 

■  E{[X  -u]  [X  Then  under  some  regularity  conditions  (see  BILLINGSLEY 

~n  ^  ~n+h 

(1968),  Theorem  20.1)  the  following  central  limit  theorem  holds  as  n  -*■ 


N (0,2)  , 

<v 


where  Z  =  Z(0)  +  2  Z(h). 

To  apply  this  result  to  simulation  output  analysis  we  need  a  method 
to  estimate  the  covariance  matrix  Z. 

<v< 

Let  f (X)  be  the  matrix-valued  spectral  density  function  for  the 
process  {X  :  n  >  0},  namely. 


f  00 


l  2(h)  e±Xh 


X  real 


Observe  that  Z  *  f (0) .  The  function  f  can  be  approximated  arbitrarily 
closely  by  the  corresponding  spectral  density  of  a  multivariate  autoregressive 
(MAR)  process  fitted  to  the  process  {X  :  n  >  0).  To  fit  a  MAR  process  of 
order  p  we  must  find  matrices  {A^  :  0  £  k  p}  so  that 

)  A,  X  .  =  Z 
,Jc  ^n-k  ^n 


where  the  Z  's  are  i.i.d.  normal  with  mean  0  and  known  convariance  Jl.  We  select 
the  order  p  of  the  MAR  process  by  applying  the  AlC-criterion  of  AKAIKE  (1976), 

(19  ).  Once  the  model  is  fitted,  it  is  a  simple  computation  to  find  the 
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spectral  density  function  of  the  MAR  process.  This  function  evaluated  at 
zero  then  provides  our  estimate  for  2. 

This  method  and  related  ones  will  be  developed  in  detail  in  JOW  (1981). 
In  Section  6  data  will  be  presented  from  an  application  of  the  method  to  the 
lake  model  presented  in  Section  3. 


6.  Application  to  the  Lake  Model 

In  this  section  we  illustrate  the  application  of  the  three  methods  to 
the  lake  model  presented  in  EXAMPLE  2.  .  Recall  the  model  is  generated  by 
the  recursion 


X  =■  PX  .  +  e  ,  n  >  0  , 

n  n-1  n  — 

where  the  e  's  are  i.i.d.  and  p  <  1.  It  can-  be  shown  that  X  =  X  as 
n  n 

n  °°  with  E{X)  *  E{Z^}  =  (l-o)  ^  E{s^}.  For  this  simulation  we  have  taken 

the  s  's  to  nave  the  distribution  P{e  <  x)  =  1  -  (l-o)e  X,  x  >  0,  which 
n  n  —  — 

results  in  X  being  exponential  with  parameter  1.  The  simulations  were 

carried  out  to  estimate  E{X).  For  the  ERM  the  return  set  A  =  [0,b] , 

where  b  =  (-log  p)/(l-p).  This  value  can  be  shown  to  maximize  the  expected 

number  of  regenerations.  The  value  p  =  0.75  was  used  which  makes  b  =  1.15. 

For  the  RBM  the  return  set  A  =  [.75,  1.25].  The  infinite  sum  of  covariance 

terms  in  the  RBM  was  truncated  at  m  (m  =  0,  5,  10,  25).  For  the  MARM  the 

2 

vector  of  observations  used  was  (X  ,  X  ,  lrrt  ,  (X  )).  A  total  of  10,000 

n  n  lu,  In  zj  n 

observations  were  generated  in  each  of  30  replications. 
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TABLE  6.1.  Simulation  Results  for  Lake  Model 
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